USE ElliDef
USE NodeInfoDef
USE metric
USE Elliptic_Solvers
IMPLICIT NONE
TYPE (NodeInfo) :: Info


INTEGER error,N
REAL (KIND=qPrec) :: xUpper(1:3),dx,dy,dz
REAL (KIND=qPrec), POINTER, DIMENSION (:,:,:) :: p,derp,u,v,w,div,pl
REAL (KIND=qPrec), ALLOCATABLE, DIMENSION (:,:) :: ut,vt,wt


INTEGER i,j,k,nx,ny,nz,l
! CHOOSE HERE THE SOLVER
!INTEGER, PARAMETER :: solver=PARDISO_SOLVER
!INTEGER, PARAMETER :: solver=IT_FOM
INTEGER, PARAMETER :: solver=DCG_SOLVER
LOGICAL iterate









! populate Info structures

ALLOCATE(Info%ElliInfo)

Info%mX=(/500,14,1/)
nx=Info%mX(1);ny=Info%mX(2);nz=Info%mX(3)
Info%nDim=2
Info%ElliInfo%ndim=Info%ndim
Info%ElliInfo%mX=Info%mX ! last dim has to be one
!!$ALLOCATE(Info%ElliInfo%p(0:Info%mX(1)+1,0:Info%mX(2)+1,0:Info%mX(3)+1))
!!$ALLOCATE(Info%ElliInfo%pold(1:Info%mX(1),1:Info%mX(2),1:Info%mX(3)))
!!$ALLOCATE(Info%ElliInfo%qe(nx+1,ny+1,nz+1,2))
!!$ALLOCATE(Info%ElliInfo%div(nx,ny,nz))

!xUpper(1:3)=(/20.d0,30.d0,1.d0/)

Info%xLower=(/10000./2200.,0.,0./)

Info%xUpper=(/15000./2200,1.,1./)
Info%dX(1:2)=(Info%xUpper(1:2)-Info%xLower(1:2))/Info%mX(1:2);Info%dX(3)=1
Info%ElliInfo%dX(1:2)=Info%dX(1:2)
Info%ElliInfo%dX(3)=Info%dX(3) ! just to have it set to something
Info%dmax=0.

!Info%mthbc=InternalBC ! assume all Neumann conditions
Info%ElliInfo%mthbc(1:6)=0
!Info%mthbc(5)=InternalBC




CALL setprob(Info)
PRINT*,'Dimension of problem',Info%mX

ALLOCATE(ut(1:ny,1:nz),vt(1:nx,1:nz))
CALL DEFINE_TOPO(Info) ! define the quantities used to calculate the metric
error=CREATE_OPERATORS(Info%ElliInfo)

! zero BC here to be safe 
DO i=1,2*Info%nDim
   Info%ElliInfo%bc(i)%p=zero
END DO

IF(.NOT.check(Info%ElliInfo,CHECK_BC)) THEN
   PRINT*,'BC are not consistent'
   STOP
END IF






! calculates divergence


u=>Info%ElliInfo%qe(1:nx+1,1:ny,1:nz,1);
v=>Info%ElliInfo%qe(1:nx,1:ny+1,1:nz,2);





CALL RANDOM_NUMBER(u)
CALL RANDOM_NUMBER(v)

CALL RANDOM_NUMBER(ut)
CALL RANDOM_NUMBER(vt)
!!$
!!$


        



! only need OP= optional to use IT wth PARDISO precon
i=POISSON_SOLVERS(Info%ElliInfo,INIT,solver)!,OP=CSR_PRECOND)
ITERATE=.true.;i=1
DO WHILE(ITERATE.AND.i<100)
PRINT*,'iteration ',i
i=i+1


Info%ElliInfo%bc(1)%p(1:ny,1:nz)=u(1,:,:)+ut ! u*-utarget
Info%ElliInfo%bc(2)%p(1:ny,1:nz)=u(nx+1,:,:)+ut
Info%ElliInfo%bc(3)%p(1:nx,1:nz)=v(:,1,:)+vt
Info%ElliInfo%bc(4)%p(1:nx,1:nz)=v(:,ny+1,:)+vt
!!$


PRINT*,'calculates divergence'


CALL take_div(Info%ElliInfo,.true.)



div=>Info%ElliInfo%div
PRINT*,'divergence stat'


!!$


PRINT*,'max divergence before = ',maxval(abs(div))

PRINT*,'rms divergence before = ',sqrt(SUM(div**2)/(nx*ny*nz))

PRINT*,'sum of div before =',SUM(div)
IF(solver==DCG_SOLVER) j=POISSON_SOLVERS(Info%ElliInfo,INIT,solver) ! need to reinit solver at each iteration
j=POISSON_SOLVERS(Info%ElliInfo,SOLVE,solver)

! subtract gradient
CALL project(Info%ElliInfo,.true.)

PRINT*,'calculates divergence'
10 CALL take_div(Info%ElliInfo,.false.)
div=>Info%ElliInfo%div

PRINT*,'max divergence after = ',maxval(abs(Info%ElliInfo%div))

PRINT*,'rms divergence after = ',sqrt(SUM(Info%ElliInfo%div**2)/(nx*ny*nz))
ITERATE=sqrt(SUM(Info%ElliInfo%div**2)/(nx*ny*nz))>1e-10

ENDDO

i=POISSON_SOLVERS(Info%ElliInfo,FINISH,solver)
PRINT*,'parison_mopup'
open(10,file='u.dat',form='unformatted')
write(10) Info%ElliInfo%qe(1:nx+1,1:ny,1:nz,1)
close(10)
open(10,file='w.dat',form='unformatted')
write(10) Info%ElliInfo%qe(1:nx,1:ny+1,1:nz,2)
close(10)
open(10,file='p.dat',form='unformatted')
write(10) Info%ElliInfo%p(1:nx,1:ny,1:nz)
close(10)


error=FREEUP_OPERATORS(Info%ElliInfo)


 
DEALLOCATE(Info%ElliInfo)
DEALLOCATE(Info%topo%h,Info%topo%hx)
DEALLOCATE(Info%topo)
DEALLOCATE(ut,vt)


END

